/* 
 * File:   Defect.h
 * Author: alrojas
 *
 * Created on 16 de junio de 2014, 10:02
 */

#ifndef DEFECT_H
#define	DEFECT_H

#include "Utils.h"
#include "cv.h"

using namespace Utils;
using namespace cv;

namespace Defect{
    
    /**
     * Erosiona una imagen usando un elemento estructurante específico
     * @param mInput Matriz de Entrada
     * @param mOutput Matriz de Salida
     * @param ksize Tamaño del kernel
     */
    inline void Eroding(Mat & src, Mat & dst, int ksize){
        Mat kernel = getStructuringElement(MORPH_RECT,Size(4 * ksize + 1, 4 * ksize + 1),Point(ksize,ksize));    
        erode(src,dst,kernel);
    }
     
    /**
     * Genera intervalo de numeros flotantes
     * @param min Valor mínimo
     * @param max Valor máximo
     * @param step Incremento
     * @return 
     */
    inline Mat CreateInterval(float min, float max, float step){

        int length = (int)(max - min + 1)/step;
        Mat interval = Mat(1,length,CV_32F);
        for(int i = 0; min <= max; i++){
            interval.at<float>(0,i) = min;
            min = min + step;
        }   
        return interval;
    }
    
    /**
     * Replica los vectores a y b para producir una malla completa
     * @param a Vector rejilla especificando una serie de puntos en el eje x
     * @param b Vector rejilla especificando una serie de puntos en el eje y
     * @param x Matriz de salida en el eje x.
     * @param y Matriz de salida en el eje y.
     */
    inline void meshgrid(Mat a, Mat b, Mat & x, Mat & y){

        x = a;
        for(int i = 1; i < b.cols; i++){
            vconcat(x,a,x);
        }

        transpose(b,b);

        y = b;
        for(int i = 1; i < a.cols; i++){
            hconcat(y,b,y);
        }
    }
    
    /**
     * Encuentra el mínimo y máximo global de una matriz.
     * @param src Entrada de un solo canal.
     * @return Punto con el valor mínimo y máximo.
     */
    inline Point MinMax_CV_32F(Mat & src){
        double min,max;
     
        minMaxLoc(src,&min,&max,0,0);

        return Point(min,max);
    }
    
    /**
     * Genera un filtro gaussiano de dos dimensiones.
     * @param dst Matriz de Salida de tamano ksize con desviación sigma.
     * @param ksize Tamano del filtro.
     * @param sigma Factor de desviación.
     */
    inline void fspecial(Mat & dst, Point ksize, float sigma){

        float e = 2.7186;
        Point size((ksize.x - 1)/2,(ksize.y - 1)/2);
        float std = sigma;

        Mat interval1 = CreateInterval(-size.x,size.x,1);
        Mat interval2 = CreateInterval(-size.y,size.y,1);
        Mat x,y;
        meshgrid(interval2,interval1,x,y);

        Mat arg, xtemp,ytemp;
        pow(x,2,xtemp);
        pow(y,2,ytemp);

        arg = -(xtemp + ytemp)/(2 * std * std);
        Mat h;
        cv::exp(arg,h);

        Point mimmax = MinMax_CV_32F(h);

        Mat htemp = (h > e * mimmax.y) ;
        h.setTo(0,htemp);

        float sumh = sum(h).val[0];

        if(sumh != 0){
            dst = h / sumh; 
         }
    }
    
    /**
     * 
     * @param src Matriz de entrada.
     * @param angle Angulo de rotación.
     * @param dst Matriz de salida.
     */
    inline void rotate(Mat src, double angle, Mat& dst){
        Mat rot_mat(2,3,CV_32F);

        Point center = Point(src.cols/2,src.rows/2);
        rot_mat = getRotationMatrix2D(center,angle,1);

        warpAffine(src,dst,rot_mat,src.size());
    }
    
    /**
     * Filtra una matriz bidimensional utilizando un filtro especifico.
     * @param src Matriz de entrada.
     * @param dst Matriz de salida.
     * @param kernel Filtro bidimensional.
     * @param conv True, para aplicar convolución. False, para no hacerlo. 
     * @param options Opciones que controlan la operación de filtro.  
     */
    inline void imfilter(Mat src,Mat & dst, Mat & kernel, bool conv, string options){

        if(conv)
            rotate(kernel,-180,kernel);

        if(options == "symmetric"){
            src.convertTo(src,CV_32F);
            filter2D(src,dst,-1,kernel,Point(-1,-1),0,BORDER_REFLECT);
        }
    }
    
    /**
     * Filtra una matriz utilizando el filtro gaussiano.
     * @param src Matriz de entrada.
     * @param dst Matriz de salida.
     * @param sigma Factor de Desviacion Estandar.
     */
    inline void imgaussian(Mat src, Mat & dst, float sigma/*, Point size*/){

        Point hsize(round(6 * sigma + 1),1);

        Mat H, Hx, Hy;
        fspecial(H,hsize,sigma);

        if(src.cols != 1 && src.rows != 1){
            transpose(H,Hx);
            Hy = H;
            Mat temp;
            imfilter(src,temp,Hy,0,"symmetric");         
            imfilter(temp,dst,Hx,0,"symmetric");
        }
        else{

        }    
    }
    
    /**
     * Replica los vectores a y b para producir una malla completa
     * @param a Vector rejilla especificando una serie de puntos en el eje x
     * @param b Vector rejilla especificando una serie de puntos en el eje y
     * @param x Matriz de salida en el eje x.
     * @param y Matriz de salida en el eje y.
     */
    inline void ndgrid(Mat a, Mat b, Mat & x, Mat & y){
    
        transpose(a,a);
        x = a;
        for(int i = 1; i < a.rows; i++){
            hconcat(x,a,x);
        }

        y = b;
        for(int i = 1; i < b.cols; i++){
            vconcat(y,b,y);
        }
    }
        
    /**
     * Filtra una matriz bidimensional utilizando un filtro especifico.
     * @param src Matriz de Entrada.
     * @param dst Matriz de Salida.
     * @param kernel Filtro bidimensional.
     */
    inline void imfilter(Mat src,Mat & dst, Mat & kernel){
        // BORDER_CONSTANT,BORDER_ISOLATED == 'conv' parameter only from MATLAB
        // BORDER_REFLECT == 'conv' + 'symmetric' parameter from MATLAB

        rotate(kernel,-180,kernel);

        src.convertTo(src,CV_32F);
        filter2D(src,dst,-1,kernel,Point(-1,-1),0,BORDER_REFLECT);
        //cout << dst << endl;
    }
    
    /**
     * Genera derivadas de la imagen basadas utilizando Gauss.
     * @param src Matriz de entrada.
     * @param dst Matriz de salida.
     * @param sigma Factor de desviación.
     * @param type Dirección.
     */
    inline void ImageDerivatives2D(Mat src, Mat & dst, float sigma, string type){

        double PI = 3.14159265; // 3.1416;
        Mat x,y;
        Mat a = CreateInterval(floor(-3 * sigma),ceil(3 * sigma),1);

        ndgrid(a,a,x,y);

        Mat m;
        if(type == "x"){ // x
            Mat xtmp = x / (2 * PI * pow(sigma,4));
            Mat ntmp,xt,yt, et;
            cv::pow(x,2,xt);
            cv::pow(y,2,yt);
            cv::add(xt,yt,ntmp);
            ntmp = -ntmp / (2 * pow(sigma,2));
            cv::exp(ntmp,et);
            cv::multiply(-xtmp,et,m);
        }
        else if(type == "y"){ // y
            Mat xtmp = y / (2 * PI * pow(sigma,4));
            Mat ntmp,xt,yt, et;
            cv::pow(x,2,xt);
            cv::pow(y,2,yt);
            cv::add(xt,yt,ntmp);
            ntmp = -ntmp / (2 * pow(sigma,2));
            cv::exp(ntmp,et);
            cv::multiply(-xtmp,et,m);
        }
        else if(type == "xx"){ // xx
            Mat ntmp,xt,yt, et;
            cv::pow(x,2,xt);
            cv::pow(y,2,yt);
            Mat xtmp = 1 / (2 * PI * pow(sigma,4)) * (xt/pow(sigma,2) - 1);
            cv::add(xt,yt,ntmp);
            ntmp = -ntmp / (2 * pow(sigma,2));
            cv::exp(ntmp,et);
            cv::multiply(xtmp,et,m);
        }
        else if(type == "xy" || type == "yx"){ // xy yx
            Mat ntmp,xt,yt, et;
            cv::pow(x,2,xt);
            cv::pow(y,2,yt);
            Mat z;
            cv::multiply(x,y,z);
            Mat xtmp = 1 / (2 * PI * pow(sigma,6)) * z;
            cv::add(xt,yt,ntmp);
            ntmp = -ntmp / (2 * pow(sigma,2));
            cv::exp(ntmp,et);
            cv::multiply(xtmp,et,m);
        }
        else if(type == "yy"){ // yy
            Mat ntmp,xt,yt, et;
            cv::pow(x,2,xt);
            cv::pow(y,2,yt);
            Mat xtmp = 1 / (2 * PI * pow(sigma,4)) * (yt/pow(sigma,2) - 1);
            cv::add(xt,yt,ntmp);
            ntmp = -ntmp / (2 * pow(sigma,2));
            cv::exp(ntmp,et);
            cv::multiply(xtmp,et,m);
            //print(m);
        }

        imfilter(src,dst,m);
    }
    
    /**
     * Calcula la energía descrita por la imagen.
     * @param src Matriz de entrada.
     * @param dst Matriz de salida.
     * @param wline Atracción a las líneas. Si es negativo para las líneas negras de otro modo blanco
     * @param wedge Atracción a los bordes.
     * @param wterm Atracción a las terminaciones; puntos finales y esquinas.
     * @param sigma Valor para calular las derivadas de la imagen.
     */
    inline void ExternalForceImage2D(Mat src, Mat & dst, float wline, float wedge, float wterm, float sigma){

        Mat x,y,xx,xy,yy;
        ImageDerivatives2D(src,x,sigma,"x");
        ImageDerivatives2D(src,y,sigma,"y");
        ImageDerivatives2D(src,xx,sigma,"xx");
        ImageDerivatives2D(src,xy,sigma,"xy");
        ImageDerivatives2D(src,yy,sigma,"yy");

        int ksize = sigma * 6;
        if((ksize % 2) == 0) ksize++;

        Mat Eline;
        imgaussian(src,Eline,sigma);

        Mat tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8;

        cv::pow(x,2.f,tmp1);
        cv::multiply(yy,tmp1,tmp2);

        cv::multiply(xy,x,tmp3);
        cv::multiply(tmp3,y,tmp4);
        tmp4 = 2 * tmp4;

        pow(y,2.f,tmp8);
        cv::multiply(xx,tmp8,tmp8);

        cv::pow(x,2.f,tmp5);
        cv::pow(y,2.f,tmp6);

        tmp7 = 1.0 + tmp5 + tmp6;

        Mat tmp20;

        cv::pow(tmp7,1.5,tmp20);

        Mat Eterm = (tmp2 - tmp4 + tmp8) / tmp7;

        Mat tmp9, tmp10;

        cv::pow(x,2.f,tmp9);
        cv::pow(y,2.f,tmp10);

        Mat Eedge;
        cv::pow(tmp9 + tmp10,0.5,Eedge);

        dst = (wline * Eline) - (wedge * Eedge) - (wterm * Eterm);    

        double minVal, maxVal;
        minMaxLoc(dst, &minVal, &maxVal);

        dst = (dst - minVal) / (maxVal - (minVal));

        dst = dst * 255; 
    }
    
    /**
     * Detecta el contorno (cuero de alpaca) de una imagen.
     * @param src Matriz de entrada.
     * @param contours Contornos detectados. Cada contorno es almacenado como un vector de puntos.
     * @param thresh Valor umbral.
     * @param ksize Tamaño del kernel.
     */
    inline void EdgeDetection(Mat & src, vector < vector < Point > > & contours, int thresh, int ksize){
        Mat grayscale = src;  
        Mat edge;
        vector < Vec4i > hierarchy;

        medianBlur(grayscale,grayscale,ksize);    
        threshold(grayscale,edge,thresh,255,THRESH_BINARY);

        floodFill(edge,Point(0,0),Scalar(0));
        floodFill(edge,Point(edge.cols - 2,0),Scalar(0));
        floodFill(edge,Point(0,edge.rows - 2),Scalar(0));
        floodFill(edge,Point(edge.cols - 2,edge.rows - 2),Scalar(0));

        findContours(edge,contours,hierarchy,CV_RETR_TREE,CV_CHAIN_APPROX_TC89_L1);
        
        uint i;
        for(i = 0; hierarchy[i][2] == -1 && hierarchy[i][3] == -1; i++);
        
        contours.erase(contours.begin(),contours.begin() + i);
          
    }
    
    /**
     * Detecta los defectos (en cueros de alpaca) de una imagen.
     * @param src Matriz de entrada.
     * @param defects Defectos detectados. El contorno es almacenado como un vector de puntos.
     * @param eKSize Tamaño kernel de erosión.
     * @param thresh Valor umbral.
     * @param edKSize Tamaño kernel de detección de contorno. 
     */
    inline void EdgeDefectsDetection(Mat & src, vector < vector < Point > > & defects, int eKSize, int thresh, int edKSize){

        vector < Vec4i > hierarchy;

        Eroding(src,src,eKSize);

        Mat grayscale;
        cvtColor(src,grayscale,CV_BGR2GRAY);  

        Mat threshold_otsu;

        GaussianBlur(grayscale,threshold_otsu,Size(edKSize,edKSize),thresh,3*thresh);
        threshold(threshold_otsu,threshold_otsu,thresh,255,THRESH_OTSU);

        findContours(threshold_otsu,defects,hierarchy,CV_RETR_TREE,CV_CHAIN_APPROX_TC89_L1,Point(0,0));

        Mat mask(src.rows,src.cols,CV_8U,Scalar(255));

        for(uint j = 0; j < defects.size(); j++){
            drawContours(mask,defects,j,Scalar(0),CV_FILLED);
        }

        grayscale.setTo(0,mask);

        Mat output;
        ExternalForceImage2D(grayscale,output,0.04,2,0.01,10);

        Mat dst;

        output.convertTo(output,CV_8U);
        equalizeHist(output,dst);
        
        EdgeDetection(dst,defects,thresh,edKSize);    
    }
    
    /**
     * Encuentra el mínimo global de un vector.
     * @param src Vector de entrada.
     * @param index Posición del minimo valor del vector.
     * @param distance Valor de la minima distancia.
     */
    inline void Min(vector < float > & src, int &index, float & distance){
        distance = src[0];
        index = 0;
         for(int i = 1; i < (int)src.size(); i++){
             if(src[i] < distance){
                 distance = src[i];
                 index = i;
             }
         }
     }
    
    /**
     * Calcula la minima distancia de un punto a un grupo de puntos.
     * @param cluster Grupo de puntos.
     * @param point Punto de referencia.
     * @return Distancia mínina.
     */
    inline float DistanceToCluster(vector < Point > & cluster, Point & point){
        float min_distance = 500;
        float d;
        for(uint i = 0; i < cluster.size(); i++){
            //d = Utils::EuclidenDistance(Point(defect.x,defect.y),Point(cluster[i].x,cluster[i].y));
            d = EuclideanDistance(point,cluster[i]);
                    
            if(d < min_distance){
                min_distance = d;
            }
        }
        return min_distance;
    }
    
    /**
     * Agrupa en cluster (grupos) un grupos de defectos.
     * @param defects Vector de defectos. 
     * @param clusters Vector de clusters (grupos).
     * @param distance Umbral de distancia.
     */
    inline void UnsupervisedClustering(vector < Point > & defects, vector < vector < Point > > & clusters, float distance){
    
        vector < Point > cluster;
        cluster.push_back(defects[0]);
        
        clusters.push_back(cluster);
        
        float d;

        for(uint i = 1; i < defects.size(); i++){
            vector < float > ds;
            for(uint j = 0; j < clusters.size(); j++){
                d = DistanceToCluster(clusters[j],defects[i]);
                
                ds.push_back(d);
            }
            
            int min_index;
            float min_distance;
            
            Min(ds,min_index,min_distance);
            
            if(min_distance <= distance){
                clusters[min_index].push_back(defects[i]);
            }
            else{
                vector < Point > new_cluster;
                new_cluster.push_back(defects[i]);
                clusters.push_back(new_cluster);
            }            
        } 
    }
    
    /**
     * Calcula el área del contorno.
     * @param contour Vector de entrada que representa el contorno.
     * @return Area del contorno.
     */
    inline float CalculateTotalArea(vector < Point > src){

        return contourArea(src,false);
    }

    /**
     * Calcula el área de la zonas defectuosas.
     * @param defects Vector de defectos.
     * @return Suma total de las zonas defectosas.
    */
    inline float CalculateDefectsArea(vector < vector < Point > > defects){
        float area;
        for(uint i = 1; i < defects.size(); i++){
            area += contourArea(defects[i],false);
        }
        return area;
    }        
} 

#endif	/* DEFECT_H */

